IE582Project

Introduction

In this project, we will do a study on data analysis and applications on 30 different stocks in Borsa Istanbul, their prices and corresponding time periods. Additionally, starting from December 24th, we will make a 10-hour stock price prediction for 30 stocks (300 in total) every day. We will analyze these predictions with different Data Mining/Machine Learning algorithms and examine which one might be more efficient to use in this report. Since there are not much raw columns in data, we need to somehow obtain an extra information from existing features or from another sources (such as Yahoo Finance using Quantmod library).

The data includes hourly average prices of selected stocks from Borsa İstanbul, presented in separate csv files for each time interval, encompassing about 4 years of data in long format. We will going to find ways to handle this long format since we do not want other stocks affects a specific stock’s price in negative way.

Let’s start with initialization.

Loading required libraries & Data.

At first, We loaded the data to the R Environment to complete the parts of the question. Also, loaded some libraries providing the visualizations to interpret the data.

library(data.table)
library(quantmod) # for Yahoo Finance library
library(ggplot2)
library(dplyr)
library(lubridate)
library(plotly)
library(rpart.plot)
library(skimr)
library(GGally)
library(caret) # for machine learning and predictive modeling
library(rpart) # for rpart
library(gbm) # for Boosted Decision Trees
library(randomForest) # for random forests
library(reshape2) # for graphs
library(FNN)
library(class)
require(rattle)
library(pROC)
library(pdp)
library(MLmetrics)
library(magick)
library(forecast)
library(cluster)
library(plotly)
library(factoextra)
library(Metrics)

Note that the path is set to the folder in my computer. Please change it to run the commands properly.

As can be seen that there are 24 different csv files including historical stock prices of 30 stocks from 2018 to 2023. With the script below, they will be merged to a single data frame.

path <- "C:/Users/anil.turgut/Desktop/582Project/Data"

# Get a list of CSV files in the directory
files <- list.files(path, pattern = "\\.csv", full.names = TRUE)
files

##  [1] "C:/Users/anil.turgut/Desktop/582Project/Data/20180101_20180401_bist30.csv"
##  [2] "C:/Users/anil.turgut/Desktop/582Project/Data/20180402_20180701_bist30.csv"
##  [3] "C:/Users/anil.turgut/Desktop/582Project/Data/20180702_20180930_bist30.csv"
##  [4] "C:/Users/anil.turgut/Desktop/582Project/Data/20181001_20181230_bist30.csv"
##  [5] "C:/Users/anil.turgut/Desktop/582Project/Data/20181231_20190331_bist30.csv"
##  [6] "C:/Users/anil.turgut/Desktop/582Project/Data/20190401_20190630_bist30.csv"
##  [7] "C:/Users/anil.turgut/Desktop/582Project/Data/20190701_20190929_bist30.csv"
##  [8] "C:/Users/anil.turgut/Desktop/582Project/Data/20190930_20191229_bist30.csv"
##  [9] "C:/Users/anil.turgut/Desktop/582Project/Data/20191230_20200329_bist30.csv"
## [10] "C:/Users/anil.turgut/Desktop/582Project/Data/20200330_20200628_bist30.csv"
## [11] "C:/Users/anil.turgut/Desktop/582Project/Data/20200629_20200927_bist30.csv"
## [12] "C:/Users/anil.turgut/Desktop/582Project/Data/20200928_20201227_bist30.csv"
## [13] "C:/Users/anil.turgut/Desktop/582Project/Data/20201228_20210328_bist30.csv"
## [14] "C:/Users/anil.turgut/Desktop/582Project/Data/20210329_20210627_bist30.csv"
## [15] "C:/Users/anil.turgut/Desktop/582Project/Data/20210628_20210926_bist30.csv"
## [16] "C:/Users/anil.turgut/Desktop/582Project/Data/20210927_20211226_bist30.csv"
## [17] "C:/Users/anil.turgut/Desktop/582Project/Data/20211227_20220327_bist30.csv"
## [18] "C:/Users/anil.turgut/Desktop/582Project/Data/20220328_20220626_bist30.csv"
## [19] "C:/Users/anil.turgut/Desktop/582Project/Data/20220627_20220925_bist30.csv"
## [20] "C:/Users/anil.turgut/Desktop/582Project/Data/20220926_20221225_bist30.csv"
## [21] "C:/Users/anil.turgut/Desktop/582Project/Data/20221226_20230326_bist30.csv"
## [22] "C:/Users/anil.turgut/Desktop/582Project/Data/20230327_20230625_bist30.csv"
## [23] "C:/Users/anil.turgut/Desktop/582Project/Data/20230626_20230924_bist30.csv"
## [24] "C:/Users/anil.turgut/Desktop/582Project/Data/20230925_20231120_bist30.csv"
## [25] "C:/Users/anil.turgut/Desktop/582Project/Data/20231121_20231223_bist30.csv"

# Create an empty data frame to store the merged data
data <- data.frame()
# Loop through each file and merge it into the 'data' data frame
for (file in files) {
  # Read the CSV file
  current_data <- read.csv(file, header = TRUE)
  
  # Merge the current data with the existing merged data
  data <- bind_rows(data, current_data)
  
}

head(data)

##                   timestamp price short_name
## 1 2018-01-02 09:00:00+03:00 15.79      THYAO
## 2 2018-01-02 10:00:00+03:00 16.01      THYAO
## 3 2018-01-02 11:00:00+03:00 16.05      THYAO
## 4 2018-01-02 12:00:00+03:00 16.05      THYAO
## 5 2018-01-02 13:00:00+03:00 16.06      THYAO
## 6 2018-01-02 14:00:00+03:00 16.05      THYAO

dim(data)

## [1] 446181      3

We have successfully loaded data as a single dataframe, now move on with the Feature Engineering part.

Feature Engineering & Yahoo Finance Library

We have three features in our dataset which are timestamp, short_name, price. Since these information are not sufficient to create strong learning models, we need to obtain new features using existing ones or from another source.

Also, there are some data modifications that will helps us in further analysis or modeling in below script.

# Changing timestamp's type to POSIX to use as time series feature
data$timestamp <- as.POSIXct(data$timestamp, format = "%Y-%m-%d %H:%M:%S", tz = "UTC")
class(data$timestamp)

## [1] "POSIXct" "POSIXt"

# Create a new feature which will be used as a key while merging further information
data$time <- as.Date(format(data$timestamp, "%Y/%m/%d"))


# Extract year, month, day, and hour
data$year <- year(data$timestamp)
data$month <- month(data$timestamp)
data$day <- day(data$timestamp)
data$hour <- hour(data$timestamp)

In below we add a new feature called “forecast_ma_6” which basically includes the Moving Average forecast with 6 records of each stock price. If there are no previous 6 records in the current record, then its forecast will be set as its original price. We except that this new feature contributes while learning algorithms

stock_symbols <- c(
  "THYAO", "AKBNK", "ARCLK", "ASELS", "BIMAS", "DOHOL", "EKGYO", "EREGL", "FROTO", "GUBRF",
  "GARAN", "KRDMD", "KCHOL", "KOZAL", "KOZAA", "PGSUS", "PETKM", "SAHOL", "SASA", "SISE",
  "TAVHL", "TKFEN", "TUPRS", "TTKOM", "TCELL", "HALKB", "ISCTR", "VAKBN", "VESTL", "YKBNK"
)

add_forecast_column <- function(data) {
  data <- data %>%
    arrange(time) %>%
    mutate(`forecast_ma_6` = ifelse(row_number() < 7, price,
                                     (lag(price, 1) + lag(price, 2) + lag(price, 3) +
                                        lag(price, 4) + lag(price, 5) + lag(price, 6)) / 6))
  return(data)
}

# Apply the function for each stock symbol
result_list <- lapply(stock_symbols, function(symbol) {
  symbol_data <- data %>% filter(short_name == symbol)
  symbol_data <- add_forecast_column(symbol_data)
  return(symbol_data)
})

# Combine the results into a single data frame
data <- bind_rows(result_list)

# Moving price column to the last place of dataframe.
data <- data[, c(setdiff(names(data), "price"), "price")]

head(data,8)

##             timestamp short_name       time year month day hour forecast_ma_6
## 1 2018-01-02 09:00:00      THYAO 2018-01-02 2018     1   2    9      15.79000
## 2 2018-01-02 10:00:00      THYAO 2018-01-02 2018     1   2   10      16.01000
## 3 2018-01-02 11:00:00      THYAO 2018-01-02 2018     1   2   11      16.05000
## 4 2018-01-02 12:00:00      THYAO 2018-01-02 2018     1   2   12      16.05000
## 5 2018-01-02 13:00:00      THYAO 2018-01-02 2018     1   2   13      16.06000
## 6 2018-01-02 14:00:00      THYAO 2018-01-02 2018     1   2   14      16.05000
## 7 2018-01-02 15:00:00      THYAO 2018-01-02 2018     1   2   15      16.00167
## 8 2018-01-02 16:00:00      THYAO 2018-01-02 2018     1   2   16      16.05500
##   price
## 1 15.79
## 2 16.01
## 3 16.05
## 4 16.05
## 5 16.06
## 6 16.05
## 7 16.11
## 8 16.13

We are going to use the quantmod library to extract related-meaningful information about our 30 stocks such as volume of the stock for a day. We will analyze whether these extracted features contribute our models or not.

Basically, quantmod is an R package that provides a framework for quantitative financial modeling and trading. It provides a rapid prototyping environment that makes modeling easier by removing the repetitive workflow issues surrounding data management and visualization.

Again, for the specified 30 stocks, start_date and end_date (must be t + 1 since it accepts as exclusive rhs), we are going to extract the information about stocks such as open, close prices, volume etc.

Also, since the quantmod library does not allow us to fetch hourly data more than 7 days, we move on with the daily average volume data of stocks.

% Note that end_date must be updated after new data has been obtained.

# Specify the start and end dates
start_date <- as.Date("2018-01-02")
end_date <- as.Date("2023-11-22")

# Create an empty data frame to store the data
all_stock_data <- data.frame()

# Loop through each stock symbol and retrieve data
for (stock_symbol in stock_symbols) {
  # Append ".IS" to the stock symbol
  full_symbol <- paste0(stock_symbol, ".IS")
  
  # Download stock data
  getSymbols(full_symbol, src = "yahoo", from = start_date, to = end_date)
  
  # Extract hourly data for the specified time range
  hourly_data <- to.hourly(get(full_symbol))
  
  # Append the data to the 'all_stock_data' data frame
  all_stock_data <- rbind(all_stock_data, data.frame(Symbol = stock_symbol, hourly_data))
}

# Remove the prefix "get.full_symbol.." from column names
cleaned_colnames <- gsub("^get.*\\.", "", colnames(all_stock_data))
colnames(all_stock_data) <- cleaned_colnames

all_stock_data$Date <- as.Date(rownames(all_stock_data))
# Print the first few rows of the combined data
head(all_stock_data)

##            Symbol  Open  High   Low Close   Volume Adjusted       Date
## 2018-01-02  THYAO 15.79 16.18 15.67 16.08 44253261    16.08 2018-01-02
## 2018-01-03  THYAO 16.05 16.23 15.92 16.20 40100544    16.20 2018-01-03
## 2018-01-04  THYAO 16.20 16.51 16.18 16.29 54106099    16.29 2018-01-04
## 2018-01-05  THYAO 16.44 16.55 16.18 16.33 55516943    16.33 2018-01-05
## 2018-01-08  THYAO 16.53 16.55 16.29 16.29 28138834    16.29 2018-01-08
## 2018-01-09  THYAO 16.29 16.45 16.10 16.10 34728669    16.10 2018-01-09

Now, we have obtained the 30 stocks’ information between start and end dates. Since we have the 10-hour information of stock prices, we can use the volume information.

To add the volume information to our previous main data, we can use merge() function since there are significant number of records, it would be extremely efficient. Let’s merge them on stock_name and the date (ex. 2022-11-10). Also, we are going to eliminate some of the features such as time from the merged final dataframe.

Moreover, we do not want missing or anomaly volume values in our dataframe. As can be seen below, there are no missing values, luckily. But, we have seen that all stock’s volumes are 0 in a specific day (The day that COVID starts in Turkey). Thus, we are going to eliminate the records whose volume is 0 to inhibit the bias.

# Merge the two dataframes based on the 'short_name' and 'Date' columns
merged_data <- merge(data, all_stock_data[, c("Symbol", "Volume", "Date")], 
                     by.x = c("short_name", "time"), by.y = c("Symbol", "Date"), all.x = TRUE)

# Rename the merged 'Volume' column to 'volume' in the 'data' dataframe
merged_data <- rename(merged_data, volume = Volume)
merged_data <- merged_data %>% arrange(short_name, timestamp)

final_data <- merged_data %>%
  select(timestamp, short_name, year, month, day, hour, volume,"forecast_ma_6", price)

head(final_data)

##             timestamp short_name year month day hour   volume forecast_ma_6
## 1 2018-01-02 09:00:00      AKBNK 2018     1   2    9 22346109        6.9475
## 2 2018-01-02 10:00:00      AKBNK 2018     1   2   10 22346109        7.0602
## 3 2018-01-02 11:00:00      AKBNK 2018     1   2   11 22346109        7.0954
## 4 2018-01-02 12:00:00      AKBNK 2018     1   2   12 22346109        7.0814
## 5 2018-01-02 13:00:00      AKBNK 2018     1   2   13 22346109        7.1024
## 6 2018-01-02 14:00:00      AKBNK 2018     1   2   14 22346109        7.1377
##    price
## 1 6.9475
## 2 7.0602
## 3 7.0954
## 4 7.0814
## 5 7.1024
## 6 7.1377

colSums(sapply(all_stock_data, is.na))

##   Symbol     Open     High      Low    Close   Volume Adjusted     Date 
##        0        0        0        0        0        0        0        0

final_data <- final_data %>%
  filter(volume != 0)

head(final_data,8)

##             timestamp short_name year month day hour   volume forecast_ma_6
## 1 2018-01-02 09:00:00      AKBNK 2018     1   2    9 22346109      6.947500
## 2 2018-01-02 10:00:00      AKBNK 2018     1   2   10 22346109      7.060200
## 3 2018-01-02 11:00:00      AKBNK 2018     1   2   11 22346109      7.095400
## 4 2018-01-02 12:00:00      AKBNK 2018     1   2   12 22346109      7.081400
## 5 2018-01-02 13:00:00      AKBNK 2018     1   2   13 22346109      7.102400
## 6 2018-01-02 14:00:00      AKBNK 2018     1   2   14 22346109      7.137700
## 7 2018-01-02 15:00:00      AKBNK 2018     1   2   15 22346109      7.070767
## 8 2018-01-02 16:00:00      AKBNK 2018     1   2   16 22346109      7.100100
##    price
## 1 6.9475
## 2 7.0602
## 3 7.0954
## 4 7.0814
## 5 7.1024
## 6 7.1377
## 7 7.1235
## 8 7.1235

colSums(sapply(final_data, is.na))

##     timestamp    short_name          year         month           day 
##             0             0             0             0             0 
##          hour        volume forecast_ma_6         price 
##             0             0             0             0

volume is successfully added to the final dataframe, also “Year,Month,Day,Hour” info are extracted from the timestamp feature. There are no missing values. We are ready to move on.

In this part, we need to create 30 different dataframe from the main dataframe since each analysis & modelling will be performed independently from another stock. Thus, we will be using these stocks separately.

thyao_data <- final_data %>% filter(short_name == "THYAO")
akbnk_data <- final_data %>% filter(short_name == "AKBNK")
arclk_data <- final_data %>% filter(short_name == "ARCLK")
asels_data <- final_data %>% filter(short_name == "ASELS")
bimas_data <- final_data %>% filter(short_name == "BIMAS")
dohol_data <- final_data %>% filter(short_name == "DOHOL")
ekgyo_data <- final_data %>% filter(short_name == "EKGYO")
eregl_data <- final_data %>% filter(short_name == "EREGL")
froto_data <- final_data %>% filter(short_name == "FROTO")
gubrf_data <- final_data %>% filter(short_name == "GUBRF")
garan_data <- final_data %>% filter(short_name == "GARAN")
krdmd_data <- final_data %>% filter(short_name == "KRDMD")
kchol_data <- final_data %>% filter(short_name == "KCHOL")
kozal_data <- final_data %>% filter(short_name == "KOZAL")
kozaa_data <- final_data %>% filter(short_name == "KOZAA")
pgsus_data <- final_data %>% filter(short_name == "PGSUS")
petkm_data <- final_data %>% filter(short_name == "PETKM")
sahol_data <- final_data %>% filter(short_name == "SAHOL")
sasa_data <- final_data %>% filter(short_name == "SASA")
sise_data <- final_data %>% filter(short_name == "SISE")
tahvl_data <- final_data %>% filter(short_name == "TAVHL")
tkfen_data <- final_data %>% filter(short_name == "TKFEN")
tuprs_data <- final_data %>% filter(short_name == "TUPRS")
ttkom_data <- final_data %>% filter(short_name == "TTKOM")
tcell_data <- final_data %>% filter(short_name == "TCELL")
halkb_data <- final_data %>% filter(short_name == "HALKB")
isctr_data <- final_data %>% filter(short_name == "ISCTR")
vakbn_data <- final_data %>% filter(short_name == "VAKBN")
vestl_data <- final_data %>% filter(short_name == "VESTL")
ykbnk_data <- final_data %>% filter(short_name == "YKBNK")

Descriptive Analysis

In this part, we are going to understand and analyze the data more by using summary, aggregated functions, sense visualizations etc.

We have data with one POSIX, one CHR and rest either int or num features, which is nice to user. And by looking at the statistical metrics, there are no obvious anomaly/extreme value in the data.

Let’s also plot some visualizations to understand or demonstrate data more. Stock Price Distribution as a boxplot and average price by months in years are also shown below.

dim(final_data)

## [1] 439007      9

str(final_data)

## 'data.frame':    439007 obs. of  9 variables:
##  $ timestamp    : POSIXct, format: "2018-01-02 09:00:00" "2018-01-02 10:00:00" ...
##  $ short_name   : chr  "AKBNK" "AKBNK" "AKBNK" "AKBNK" ...
##  $ year         : num  2018 2018 2018 2018 2018 ...
##  $ month        : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ day          : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ hour         : int  9 10 11 12 13 14 15 16 17 18 ...
##  $ volume       : num  22346109 22346109 22346109 22346109 22346109 ...
##  $ forecast_ma_6: num  6.95 7.06 7.1 7.08 7.1 ...
##  $ price        : num  6.95 7.06 7.1 7.08 7.1 ...

summary(final_data)

##    timestamp                       short_name             year     
##  Min.   :2018-01-02 09:00:00.00   Length:439007      Min.   :2018  
##  1st Qu.:2019-06-21 12:00:00.00   Class :character   1st Qu.:2019  
##  Median :2020-12-10 17:00:00.00   Mode  :character   Median :2020  
##  Mean   :2020-12-11 03:22:19.56                      Mean   :2020  
##  3rd Qu.:2022-06-01 16:00:00.00                      3rd Qu.:2022  
##  Max.   :2023-11-21 18:00:00.00                      Max.   :2023  
##      month             day             hour           volume         
##  Min.   : 1.000   Min.   : 1.00   Min.   : 9.00   Min.   :9.000e+00  
##  1st Qu.: 3.000   1st Qu.: 8.00   1st Qu.:11.00   1st Qu.:7.063e+06  
##  Median : 6.000   Median :16.00   Median :13.00   Median :3.042e+07  
##  Mean   : 6.452   Mean   :15.63   Mean   :13.49   Mean   :6.474e+07  
##  3rd Qu.: 9.000   3rd Qu.:23.00   3rd Qu.:16.00   3rd Qu.:8.191e+07  
##  Max.   :12.000   Max.   :31.00   Max.   :18.00   Max.   :1.073e+09  
##  forecast_ma_6          price         
##  Min.   :  0.5844   Min.   :  0.5817  
##  1st Qu.:  5.0217   1st Qu.:  5.0229  
##  Median : 11.2297   Median : 11.2318  
##  Mean   : 32.0485   Mean   : 32.0766  
##  3rd Qu.: 26.4633   3rd Qu.: 26.4800  
##  Max.   :961.7500   Max.   :975.0000

ggplot(final_data, aes(x = short_name, y = price, fill = short_name)) +
  geom_boxplot() +
  labs(title = "Stock Price Distribution",
       x = "Stock",
       y = "Price") +
  theme_minimal() +
  theme(legend.position = "none")

average_price_by_month <- aggregate(price ~ month + year, data = final_data, mean)

ggplot(average_price_by_month, aes(x = month, y = price, fill = factor(year))) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Average Stock Price by Month",
       x = "Month",
       y = "Average Price") +
  theme_minimal() +
  theme(legend.position = "top")

We have write a function detecting anomaly in the plots. Anomaly is detected in date, if the value is less or grater than the mean (data) +- 3 * stdev(data). There are 5 example stock prices plots to detect anomalies. None of them included any anomaly values (Rest of them are also controlled). This is also a good sign for us to create sense models.

plot_anomaly <- function(df){
  anomaly_threshold_price <-  mean(df$price) - 3 * sqrt(var(df$price)) # Adjust this threshold as needed
  anomaly_values_price <- df$price < anomaly_threshold_price
  
  stock_title <- paste(df$short_name[1],"Stock Prices Over Time with Anomalies", sep = " ")
  
  ggplot(df, aes(x = timestamp, y = price)) +
    geom_line() +
    geom_point(data = df[anomaly_values_price, ], aes(color = "Anomaly"), size = 3) +
    labs(x = "Datetime", y = "Stock Closing Prices", title = stock_title) +
    scale_color_manual(values = c("Anomaly" = "red")) +
    theme_minimal()
}
plot_anomaly(thyao_data)

plot_anomaly(akbnk_data)

plot_anomaly(vestl_data)

plot_anomaly(froto_data)

plot_anomaly(sise_data)

Finally, using quantmod library and its functionalities lets look at the visualizations of example stocks information such as THYAO.IS and SISE.IS. These candlestick and line plots are helpful for further analysis also since we get insight how current volume affects further stock price etc.

Also, all stock prices are plotted as single plot in below. It can be seen that some of the stocks have multiplied prices accross many and we think that the splitting dataset as different stocks might be good idea to block any negative effects.

# Plot candlestick chart for one stock (e.g., THYAO, SISE)
thYaoData <- getSymbols("THYAO.IS", src = "yahoo", from = "2018-01-02", to = "2023-11-22", auto.assign = FALSE)
candleChart(thYaoData, theme = "white")

siseData <- getSymbols("SISE.IS", src = "yahoo", from = "2018-01-02", to = "2023-11-22", auto.assign = FALSE)
candleChart(siseData, theme = "white")

ggplot(all_stock_data, aes(x = Date, y = Close, color = Symbol)) +
  geom_line() +
  labs(title = "Closing Prices Over Time for Selected Stocks",
       x = "Date",
       y = "Closing Price",
       color = "Stock Symbol") +
  theme_minimal()

Let’s also analyze the data by applying clustering methods. It will help us to see which stocks are similar to each other in terms of distance.

Clustering Analysis

Clustering is the task of dividing the population or data points into a number of groups such that data points in the same groups are more similar to other data points in the same group and dissimilar to the data points in other groups. It is basically a collection of objects on the basis of similarity and dissimilarity between them.

In this dataset, we can obtain insights about which stocks are similar etc.

Let’s start with k means clustering :

cluster_data <- final_data

features <- cluster_data[, c("volume", "forecast_ma_6", "price")]

# Scale the features
scaled_features <- scale(features)

# Choose the number of clusters (k)
k <- 4

# Apply k-means clustering
kmeans_result <- kmeans(scaled_features, centers = k, nstart = 25)

# Add the cluster labels to your data
cluster_data$cluster <- as.factor(kmeans_result$cluster)

head(cluster_data)

##             timestamp short_name year month day hour   volume forecast_ma_6
## 1 2018-01-02 09:00:00      AKBNK 2018     1   2    9 22346109        6.9475
## 2 2018-01-02 10:00:00      AKBNK 2018     1   2   10 22346109        7.0602
## 3 2018-01-02 11:00:00      AKBNK 2018     1   2   11 22346109        7.0954
## 4 2018-01-02 12:00:00      AKBNK 2018     1   2   12 22346109        7.0814
## 5 2018-01-02 13:00:00      AKBNK 2018     1   2   13 22346109        7.1024
## 6 2018-01-02 14:00:00      AKBNK 2018     1   2   14 22346109        7.1377
##    price cluster
## 1 6.9475       4
## 2 7.0602       4
## 3 7.0954       4
## 4 7.0814       4
## 5 7.1024       4
## 6 7.1377       4

# Scatter plot in 3D
plot_ly(cluster_data, x = ~volume, y = ~`forecast_ma_6`, z = ~price, color = ~cluster,
        text = ~short_name,  # This adds stock names as labels
        type = "scatter3d", mode = "markers", marker = list(size = 5)) %>%
  layout(scene = list(xaxis = list(title = "Volume"),
                      yaxis = list(title = "forecast_ma_6"),
                      zaxis = list(title = "Price"),
                      margin = list(l = 0, r = 0, b = 0, t = 0)))

fviz_cluster(kmeans_result, data = cluster_data[, c("volume", "price")])

set.seed(10)  # Set seed for reproducibility
sample_size <- 40  # Adjust the sample size as needed
subset_data <- cluster_data[sample(1:nrow(cluster_data), sample_size), 
                            c("short_name","volume", "forecast_ma_6", "price")]
normalized_data <- scale(subset_data[, -1])

distance_matrix <- dist(normalized_data, method = "euclidean")
hierarchical_cluster <- hclust(distance_matrix, method = "average")
clusters <- cutree(hierarchical_cluster, k = 5)  # You can adjust the number of clusters (k) as needed,
dend <- as.dendrogram(hierarchical_cluster)

# Plot the hierarchical clustering with factoextra
fviz_dend(dend, k = 5, main = "Hierarchical Clustering Dendrogram")

It can be seen that some stocks are stored in same clusters with minimum distances:

• VKBANK and HLKBNK

• THYAO, SISE, SASA, EREGL

• ASELS, TUPRS, YKBANK, AKBANK

These stocks can be used to predict the price trend of each other.

We are ready to move on with the modelling part of the report.

Machine Learning / Data Mining Models

We have completed the preprocessing, feature engineering, cleaning and the overall descriptive analysis part of the work. Now, we are going to use different kind of machine learning and statistical models to predict the future forecasts well.

We plan to implement 4 main models in our work :

• Decision Tree Algorithm

• Random Forest Algorithm

• Gradient Boosting Machines (GBM)

• ARIMA (AutoRegressive Integrated Moving Average)

Note that our main metric will be Weighted Mean Absolute Percentage Error (WMAPE) which is a measure used to evaluate the accuracy of forecast models. This metric provides a normalized error rate across all predictions, allowing for a fair comparison between different forecasting models or different data sets.

In this report, we will be analyzing a stock to select the best model for predictions. But in our work, we test with multiple stocks to get balanced information. We are going to use “THYAO” stock as a sample stock.

Let’s start with the Decision Tree Algorithm.

Decision Tree Algorithm

In order not to mutate the original data, we copy it and split it into training and test dataset which will be used in modeling and testing respectively.

dt_thyao_data <- thyao_data

validationIndex <- createDataPartition(dt_thyao_data$price, p=0.70, list=FALSE)

dt_train <- dt_thyao_data[validationIndex,] # 70% of data to training
dt_test <- dt_thyao_data[-validationIndex,]

cat("Dimension of the main dataset:",dim(dt_thyao_data))

## Dimension of the main dataset: 14636 9

cat("Dimension of the train dataset:",dim(dt_train))

## Dimension of the train dataset: 10247 9

cat("Dimension of the test dataset:",dim(dt_test))

## Dimension of the test dataset: 4389 9

Dataset is ready to be learned and now we are going to implement a manual cross validation to tune the hyperparameters for the decision tree and we are going to look at its MSE value in each iteration.

Assuming the tuning hyperparameter will be only the minimal number of observations per tree leaf. Setting complexity parameter to zero and minimum number of observations to split as the twice as the minimal number of observations per tree leaf.

In the end, we are going to store them in df to find the best model.

cv_results <- data.frame(minbucket = numeric(), RMSE = numeric())

set.seed(10)

# Define a range of minsplit values to try
minbucket_values <- c(1, 5, 10, 15, 20) 

# Perform cross-validation for each minbucket value
for (minbucket_val in minbucket_values) {
  # Create a decision tree regression model with the current minbucket value
  tree_model <- rpart(price~.,dt_train,method='anova',
                      control=rpart.control(cp=0, minbucket = minbucket_val, minsplit = 2*minbucket_val))
  
  predictions <- predict(tree_model, newdata = dt_train)
  
  true_values <- dt_train$price
  
  rmse <- round(sqrt(mean((predictions - true_values)^2)),3)
  
  # Store the results in the data frame
  cv_results <- rbind(cv_results, data.frame(minbucket = minbucket_val, RMSE = rmse))
}

ggplot(cv_results, aes(x = minbucket, y = RMSE)) +
  geom_line() +
  geom_point() +
  labs(title = "RMSE vs. MinBucket",
       x = "minbucket",
       y = "RMSE") +
  theme_minimal()

Least MSE belongs to the minbucket = 1 but the this model poorly performs in the test dataset. It seems that overfitting occurs. On the other hand, minbucket = 5 performs well both in training and test dataset. Looks like the model with minbucket = 5 is the best model in terms of RMSE.

In below, we have created the best model with corresponding hyperparameters and some visualizations are plotted.

minbucket_val <- 5
best_tree_model <- rpart(price~.,dt_train,method='anova',
                         control=rpart.control(cp=0, minbucket = minbucket_val, minsplit = 2*minbucket_val))

fancyRpartPlot(best_tree_model)

barplot(best_tree_model$variable.importance,horiz=F)

Lets look at its other regression metrics such as MAE, RMSE, MAPE etc. Also predicted vs. actual plot can be seen:

# Make predictions on the test set
predictions <- predict(best_tree_model, newdata = dt_test)
# True values from the test set
true_values <- dt_test$price

# Calculate Mean Absolute Error (MAE)
mae <- mean(abs(predictions - true_values))
print(paste("Mean Absolute Error (MAE):", mae))

## [1] "Mean Absolute Error (MAE): 0.46083640197174"

# Calculate Root Mean Squared Error (RMSE)
rmse <- round(sqrt(mean((predictions - true_values)^2)),3)
print(paste("Root Mean Squared Error (RMSE):", rmse))

## [1] "Root Mean Squared Error (RMSE): 1.114"

# Calculate Mean Absolute Percentage Error (MAPE)
mape <- mean(abs((predictions - true_values) / true_values) * 100, na.rm = TRUE)
cat("Mean Absolute Percentage Error (MAPE):", mape, "\n")

## Mean Absolute Percentage Error (MAPE): 0.9371201

calculate_wmape <- function(actual, forecasted) {
  n <- length(actual)
  wmape <- (sum(abs(actual - forecasted)) / sum(abs(actual))) * 100
  return(wmape)
}

# Calculate Weighted Mean Absolute Percentage Error (WMAPE)
wmape <- calculate_wmape(predictions, true_values)
print(paste("Weighted Mean Absolute Percentage Error (WMAPE):", wmape))

## [1] "Weighted Mean Absolute Percentage Error (WMAPE): 0.969655369325999"

# Calculate R-squared (R²)
sse <- sum((predictions - true_values)^2)
sst <- sum((true_values - mean(true_values))^2)
rsquared <- 1 - (sse / sst)
print(paste("R-squared (R²):", rsquared))

## [1] "R-squared (R²): 0.999701935778422"

# Create a scatterplot
plot(true_values, predictions, 
     main = "Predictions vs Actual Data",
     xlab = "True Values",
     ylab = "Predicted Values",
     pch = 16,  # Set the point character
     col = "blue"  # Set the point color
)

# Add a diagonal line for reference
abline(0, 1, col = "red", lty = 2)

# Add a legend
legend("topright", legend = "Diagonal Line (Reference)", col = "red", lty = 2, cex = 0.8)

It is extremely successful with the test dataset also, but it might get overfit by the complex model since Decision Trees tend to overfit easily.

Let’s move on with Random Forest.

Random Forest Algorithm

In order not to mutate the original data, we copy it and split it into training and test dataset which will be used in modeling and testing respectively.

rf_thyao_data <- thyao_data
validationIndex <- createDataPartition(rf_thyao_data$price, p=0.70, list=FALSE)

rf_train <- rf_thyao_data[validationIndex,] # 70% of data to training
rf_test <- rf_thyao_data[-validationIndex,]

cat("Dimension of the main dataset:",dim(rf_thyao_data))

## Dimension of the main dataset: 14636 9

cat("Dimension of the train dataset:",dim(rf_train))

## Dimension of the train dataset: 10247 9

cat("Dimension of the test dataset:",dim(rf_test))

## Dimension of the test dataset: 4389 9

Dataset is ready to be learned and now we are going to implement a manual cross validation to tune the hyperparameters for the random forest and we are going to look at its error value in each iteration.

Assuming the tuning hyperparameter will be only the effect of the ratio of the number of features evaluated at each split (mtry) and setting other parameters as J=500 and the minimal number of observations per tree leaf=5.

In the end, we are going to store them in table to find the best model.

set.seed(10)
# Set the number of trees and nodesize
num_trees <- 500
min_obs_per_leaf <- 5

# Create a grid of mtry values to explore
mtry_values <- c(2, 4, 6, 8, 10)  # Add more values as needed

# Create an empty data frame to store results
rf_results <- data.frame(mtry = numeric(), RMSE = numeric())

# Perform grid search
for (m in mtry_values) {
  # Train the Random Forest model
  rf_model <- randomForest(price ~ ., data = rf_train, 
                           ntree = num_trees, nodesize = min_obs_per_leaf, mtry = m)
  
  # Make predictions on the training set
  predictions <- predict(rf_model, newdata = rf_train)
  
  true_values <- rf_train$price
  
  rmse <- round(sqrt(mean((predictions - true_values)^2)),3)
  
  # Store the results in the data frame
  rf_results <- rbind(rf_results, data.frame(mtry = m, RMSE = rmse))
}

# Print the results
print(rf_results)

##   mtry  RMSE
## 1    2 1.013
## 2    4 0.373
## 3    6 0.377
## 4    8 0.383
## 5   10 0.384

ggplot(rf_results, aes(x = mtry, y = RMSE)) +
  geom_line() +
  geom_point() +
  labs(title = "RMSE vs. mtry",
       x = "mtry",
       y = "RMSE") +
  theme_minimal()

Looks like the model converges when mtry = 4 (given ntree = 500 & nodesize = 5) is the best model in terms of RMSE, lets also analyze the other metrics using test dataset.

In below, we have created the best model with corresponding hyperparameters and metrics with test data are listed.

set.seed(10)

best_rf_model <- randomForest(price ~ ., data = rf_train, 
                              ntree = num_trees, nodesize = min_obs_per_leaf, mtry = 4)

# Make predictions on the test set
predictions <- predict(best_rf_model, newdata = rf_test)

# True values from the test set
true_values <- rf_test$price

# Calculate Mean Absolute Error (MAE)
mae <- mean(abs(predictions - true_values))
print(paste("Mean Absolute Error (MAE):", mae))

## [1] "Mean Absolute Error (MAE): 0.345629370976017"

# Calculate Root Mean Squared Error (RMSE)
rmse <- round(sqrt(mean((predictions - true_values)^2)),3)
print(paste("Root Mean Squared Error (RMSE):", rmse))

## [1] "Root Mean Squared Error (RMSE): 0.861"

# Calculate Mean Absolute Percentage Error (MAPE)
mape <- mean(abs((predictions - true_values) / true_values) * 100, na.rm = TRUE)
cat("Mean Absolute Percentage Error (MAPE):", mape, "\n")

## Mean Absolute Percentage Error (MAPE): 0.6648455

# Calculate Weighted Mean Absolute Percentage Error (WMAPE)
wmape <- calculate_wmape(predictions, true_values)
print(paste("Weighted Mean Absolute Percentage Error (WMAPE):", wmape))

## [1] "Weighted Mean Absolute Percentage Error (WMAPE): 0.71896447367975"

# Calculate R-squared (R²)
sse <- sum((predictions - true_values)^2)
sst <- sum((true_values - mean(true_values))^2)
rsquared <- 1 - (sse / sst)
print(paste("R-squared (R²):", rsquared))

## [1] "R-squared (R²): 0.999823825616835"

It can be seen that random forest performed also very well on this dataset. Moreover, the possibility for overfitting the random forest model is much less than the decision trees since it is using bootstrap aggregating technique.

Let’s also make some visualizations on the grey-box random forest model.

# Plot variable importance
varImpPlot(best_rf_model)

# Create partial dependence plot
partial_plot <- partial(best_rf_model, pred.var = 'year', data = dt_test)

# Plot the partial dependence plot
plot(partial_plot)

Let’s move on with the third algorithm which is ARIMA.

ARIMA (AutoRegressive Integrated Moving Average)

In order not to mutate the original data, we copy it and split it into training and test dataset which will be used in modeling and testing respectively.

arima_thyao_data <- thyao_data
validationIndex <- createDataPartition(arima_thyao_data$price, p=0.70, list=FALSE)

arima_train <- arima_thyao_data[validationIndex,] # 70% of data to training
arima_test <- arima_thyao_data[-validationIndex,]

cat("Dimension of the main dataset:",dim(arima_thyao_data))

## Dimension of the main dataset: 14636 9

cat("Dimension of the train dataset:",dim(arima_train))

## Dimension of the train dataset: 10247 9

cat("Dimension of the test dataset:",dim(arima_test))

## Dimension of the test dataset: 4389 9

Let’s prepare the timeseriers data as 10-hours frequency, fit the ARIMA model using auto.arima function and plot.

# frequency is 10 since we have 10-hours data for each stock and day.
ts_data <- ts(arima_train$price, frequency = 10)
# plotting time-series data with prices
plot(ts_data, main = "Time Series Plot")

diff_ts_data <- diff(ts_data)
# fitting arima model
arima_model <- auto.arima(ts_data)

checkresiduals(arima_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(5,2,0)
## Q* = 551.83, df = 15, p-value < 2.2e-16
## 
## Model df: 5.   Total lags used: 20

summary(arima_model)

## Series: ts_data 
## ARIMA(5,2,0) 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ar5
##       -0.7705  -0.6258  -0.4429  -0.3222  -0.1608
## s.e.   0.0098   0.0120   0.0127   0.0120   0.0098
## 
## sigma^2 = 0.7261:  log likelihood = -12895.75
## AIC=25803.49   AICc=25803.5   BIC=25846.9
## 
## Training set error measures:
##                        ME      RMSE      MAE          MPE      MAPE      MASE
## Training set 0.0008390498 0.8518491 0.353875 0.0006634895 0.6998595 0.2934095
##                     ACF1
## Training set -0.01490018

Let’s now forecast further along with number of test records.

forecast_values <- forecast(arima_model, h = nrow(arima_test))  
plot(forecast_values, main = "ARIMA Forecast")

combined_data <- data.frame(
  timestamp = arima_test$timestamp,
  actual = arima_test$price,
  forecast_mean = forecast_values$mean,
  forecast_lower_80 = forecast_values$lower[, "80%"],
  forecast_upper_80 = forecast_values$upper[, "80%"],
  forecast_lower_95 = forecast_values$lower[, "95%"],
  forecast_upper_95 = forecast_values$upper[, "95%"]
)

ggplot(combined_data, aes(x = timestamp)) +
  geom_line(aes(y = actual, color = "Actual"), linetype = "solid", size = 1, alpha = 0.8) +
  geom_line(aes(y = forecast_mean, color = "Forecast"), linetype = "solid", size = 1, alpha = 0.8) +
  geom_ribbon(aes(ymin = forecast_lower_80, ymax = forecast_upper_80, fill = "80% CI"), alpha = 0.3) +
  geom_ribbon(aes(ymin = forecast_lower_95, ymax = forecast_upper_95, fill = "95% CI"), alpha = 0.3) +
  labs(title = "ARIMA Forecast vs Actual",
       x = "Timestamp",
       y = "Price") +
  scale_color_manual(values = c("Actual" = "blue", "Forecast" = "red")) +
  scale_fill_manual(values = c("80% CI" = "pink", "95% CI" = "orange")) +
  theme_minimal() +
  theme(legend.position = "top")

Finally, let’s look at the metrics of ARIMA model:

mae <- mean(abs(forecast_values$mean - arima_test$price))
mse <- mean((forecast_values$mean - arima_test$price)^2)
rmse <- sqrt(mse)
mape <- mean(abs((forecast_values$mean - arima_test$price) / arima_test$price) * 100, na.rm = TRUE)
wmape <- calculate_wmape(forecast_values$mean, arima_test$price)
cat("MAE:", mae, "\n")

## MAE: 6234.672

cat("MSE:", mse, "\n")

## MSE: 50610983

cat("RMSE:", rmse, "\n")

## RMSE: 7114.14

cat("MAPE:", mape, "\n")

## MAPE: 27133.35

cat("WMAPE:", wmape, "\n")

## WMAPE: 99.26006

Results are not as good as Decision Tree or Random Forest but it does not create complex models even they are simple. So, it might be more beneficial with unseen data when it comes.

Finally, let’s also model with Gradient Boosting Machine Algorithm (GBM).

GBM Algorithm

In order not to mutate the original data, we copy it and split it into training and test dataset which will be used in modeling and testing respectively. Since GBM could not handle POSIX and single type factor data types, we are removing them from the data. It might result with the lack of info but also less complex model. Let’s analyze it.

thyao_data_gbm <- thyao_data

thyao_data_gbm <- thyao_data_gbm[, !colnames(thyao_data_gbm) %in% c("short_name","timestamp")]

validationIndex <- createDataPartition(thyao_data_gbm$price, p=0.70, list=FALSE)

gbm_train <- thyao_data_gbm[validationIndex,] # 70% of data to training
gbm_test <- thyao_data_gbm[-validationIndex,]

cat("Dimension of the main dataset:",dim(thyao_data_gbm))

## Dimension of the main dataset: 14636 7

cat("Dimension of the train dataset:",dim(gbm_train))

## Dimension of the train dataset: 10247 7

cat("Dimension of the test dataset:",dim(gbm_test))

## Dimension of the test dataset: 4389 7

In above, there are 2 train & test pairs occur. One will be used in the cross validation and the other will be used to test the best model with tuned parameters.

Let’s tune the hyperparameters of the GBM with respect to the “RMSE” value.

We are mainly interested in tuning the depth (interaction.depth), the number of trees (n.trees) and the learning rate (shrinkage).

set.seed(10)

n_folds=10

fitControl=trainControl(method = "cv",
                        number = n_folds)
## gradient boosting
gbmGrid=expand.grid(interaction.depth = c(3, 5), 
                    n.trees = c(1:5)*100, 
                    shrinkage = c(0.05,0.1),
                    n.minobsinnode = 10)
set.seed(1)                        
gbm_fit=train(price ~ ., data = gbm_train, 
              method = "gbm", 
              trControl = fitControl, metric='RMSE',
              tuneGrid = gbmGrid,
              verbose=F) #verbose is an argument from gbm, prints to screen
gbm_fit

## Stochastic Gradient Boosting 
## 
## 10247 samples
##     6 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 9223, 9223, 9223, 9221, 9221, 9222, ... 
## Resampling results across tuning parameters:
## 
##   shrinkage  interaction.depth  n.trees  RMSE      Rsquared   MAE      
##   0.05       3                  100      1.344977  0.9996400  0.7298514
##   0.05       3                  200      1.089747  0.9997121  0.5020459
##   0.05       3                  300      1.070201  0.9997226  0.4923995
##   0.05       3                  400      1.060193  0.9997279  0.4887095
##   0.05       3                  500      1.051979  0.9997324  0.4858377
##   0.05       5                  100      1.201309  0.9996999  0.6235743
##   0.05       5                  200      1.037298  0.9997397  0.4712037
##   0.05       5                  300      1.019495  0.9997487  0.4634243
##   0.05       5                  400      1.008073  0.9997544  0.4585960
##   0.05       5                  500      1.001241  0.9997577  0.4558954
##   0.10       3                  100      1.185168  0.9996604  0.5651060
##   0.10       3                  200      1.132711  0.9996896  0.5399961
##   0.10       3                  300      1.107273  0.9997034  0.5283012
##   0.10       3                  400      1.097696  0.9997087  0.5226209
##   0.10       3                  500      1.087318  0.9997140  0.5179894
##   0.10       5                  100      1.101090  0.9997070  0.5160645
##   0.10       5                  200      1.065757  0.9997258  0.5017611
##   0.10       5                  300      1.044126  0.9997367  0.4924592
##   0.10       5                  400      1.029761  0.9997439  0.4857518
##   0.10       5                  500      1.020569  0.9997484  0.4817506
## 
## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were n.trees = 500, interaction.depth =
##  5, shrinkage = 0.05 and n.minobsinnode = 10.

plot(gbm_fit)

Selecting the best model is below:

Looks like the best GBM model with the parameters (Shrinkage = 0.1, depth = 5, n.trees = 500).

Let’s analyze this model with test data finally.

noftrees=500
depth=5
learning_rate=0.1

boosting_model=gbm(price~., data=gbm_train,distribution = 'gaussian', n.trees = noftrees,
                   interaction.depth = depth, n.minobsinnode = 10, shrinkage =learning_rate, cv.folds = 10)

gbm.perf(boosting_model, method = "cv")

## [1] 497

predictions <- predict(boosting_model, newdata = gbm_test, type = "response")

# True values from the test set
true_values <- gbm_test$price

# Calculate Mean Absolute Error (MAE)
mae <- mean(abs(predictions - true_values))
print(paste("Mean Absolute Error (MAE):", mae))

## [1] "Mean Absolute Error (MAE): 0.5023089255232"

# Calculate Root Mean Squared Error (RMSE)
rmse <- round(sqrt(mean((predictions - true_values)^2)),3)
print(paste("Root Mean Squared Error (RMSE):", rmse))

## [1] "Root Mean Squared Error (RMSE): 1.116"

# Calculate Mean Absolute Percentage Error (MAPE)
mape <- mean(abs((predictions - true_values) / true_values) * 100, na.rm = TRUE)
cat("Mean Absolute Percentage Error (MAPE):", mape, "\n")

## Mean Absolute Percentage Error (MAPE): 1.142101

# Calculate Weighted Mean Absolute Percentage Error (WMAPE)
wmape <- calculate_wmape(predictions, true_values)
print(paste("Weighted Mean Absolute Percentage Error (WMAPE):", wmape))

## [1] "Weighted Mean Absolute Percentage Error (WMAPE): 1.05087377811717"

It has better results than ARIMA but worse than Decision Tree and Random Forest. But still the results are satisfied eventhough the problem is complex and the information we have is limited.

Anyway we need unseen data to select the best model for predictions.

Predictions with Selected Model for All Stocks

Assuming that GBM is the selected model and lets do the predictions on the created new data

# List of stock names
stock_names <- c("THYAO", "AKBNK", "ARCLK", "ASELS", "BIMAS", "DOHOL", "EKGYO", "EREGL", "FROTO", "GUBRF",
                 "GARAN", "KRDMD", "KCHOL", "KOZAL", "KOZAA", "PGSUS", "PETKM", "SAHOL", "SASA", "SISE",
                 "TAVHL", "TKFEN", "TUPRS", "TTKOM", "TCELL", "HALKB", "ISCTR", "VAKBN", "VESTL", "YKBNK")

# Function to generate new data for each stock
generate_new_data <- function(stock_name) {
  last_timestamp <- max(final_data$timestamp[final_data$short_name == stock_name])
  last_datetime <- as.POSIXct(last_timestamp)
  
  # Increase last_datetime by one day
  last_datetime <- last_datetime + days(1)
  
  # Extract date part
  last_date <- format(last_datetime, "%Y-%m-%d")
  
  # Combine with 09:00
  start_datetime <- as.POSIXct(paste(last_date, "09:00:00"))
  
  # Generate new timestamps starting from the next day of the last record, between 9:00 and 18:00
  new_timestamps <- seq.POSIXt(from = start_datetime, 
                               by = "1 hour", length.out = 10)
  
  # Extract year, month, day, and hour
  new_data <- data.frame(
    timestamp = new_timestamps,
    short_name = rep(stock_name, 10),
    year = format(new_timestamps, "%Y"),
    month = format(new_timestamps, "%m"),
    day = format(new_timestamps, "%d"),
    hour = format(new_timestamps, "%H"),
    volume = rep(71142665, 10),
    forecast_ma_6 = rep(NA, 10),
    price = rep(NA, 10)
  )
  
  return(new_data)
}

# List to store new data for each stock
new_data_list <- list()

# Generate new data for each stock
for (stock_name in stock_names) {
  new_data <- generate_new_data(stock_name)
  new_data_list[[length(new_data_list) + 1]] <- new_data
}

# Combine new data for all stocks into a single dataframe
new_data_df <- do.call(rbind, new_data_list)

head(new_data_df,10)

##              timestamp short_name year month day hour   volume forecast_ma_6
## 1  2023-11-22 09:00:00      THYAO 2023    11  22   09 71142665            NA
## 2  2023-11-22 10:00:00      THYAO 2023    11  22   10 71142665            NA
## 3  2023-11-22 11:00:00      THYAO 2023    11  22   11 71142665            NA
## 4  2023-11-22 12:00:00      THYAO 2023    11  22   12 71142665            NA
## 5  2023-11-22 13:00:00      THYAO 2023    11  22   13 71142665            NA
## 6  2023-11-22 14:00:00      THYAO 2023    11  22   14 71142665            NA
## 7  2023-11-22 15:00:00      THYAO 2023    11  22   15 71142665            NA
## 8  2023-11-22 16:00:00      THYAO 2023    11  22   16 71142665            NA
## 9  2023-11-22 17:00:00      THYAO 2023    11  22   17 71142665            NA
## 10 2023-11-22 18:00:00      THYAO 2023    11  22   18 71142665            NA
##    price
## 1     NA
## 2     NA
## 3     NA
## 4     NA
## 5     NA
## 6     NA
## 7     NA
## 8     NA
## 9     NA
## 10    NA

Filling the forecast_ma_6 of new_data column with using the previous data for each stock.

# Function to fill forecast_ma_6 column for each stock in new_data_df
fill_forecast_ma_6 <- function(stock_name, new_data_df, original_data) {
  # Filter original data for the specific stock
  stock_data <- original_data[original_data$short_name == stock_name, ]
  
  # Set initial values for the moving average calculation
  window_size <- 6
  forecast_column <- "forecast_ma_6"
  
  # Iterate over rows in new_data_df for the current stock
  for (i in which(new_data_df$short_name == stock_name)) {
    # Extract relevant information from new_data_df
    timestamp <- new_data_df$timestamp[i]
    
    # Filter previous records from stock_data for the current stock
    previous_records <- stock_data[stock_data$timestamp < timestamp, ]
    
    # Calculate moving average forecast
    if (nrow(previous_records) >= window_size) {
      # Use mean to calculate the average iteratively
      new_data_df[i, forecast_column] <- mean(tail(previous_records[, forecast_column], window_size))
      
      # Append a new record to stock_data
      new_record <- tail(previous_records, 1)
      new_record$timestamp <- timestamp
      new_record[, forecast_column] <- new_data_df[i, forecast_column]
      stock_data <- rbind(stock_data, new_record)
    }
  }
  
  return(new_data_df)
}
# Iterate over each stock to fill forecast_ma_6 in new_data_df
for (stock_name in stock_names) {
  new_data_df <- fill_forecast_ma_6(stock_name, new_data_df, final_data)
}

head(new_data_df,10)

##              timestamp short_name year month day hour   volume forecast_ma_6
## 1  2023-11-22 09:00:00      THYAO 2023    11  22   09 71142665      258.7014
## 2  2023-11-22 10:00:00      THYAO 2023    11  22   10 71142665      259.1030
## 3  2023-11-22 11:00:00      THYAO 2023    11  22   11 71142665      259.5021
## 4  2023-11-22 12:00:00      THYAO 2023    11  22   12 71142665      259.7872
## 5  2023-11-22 13:00:00      THYAO 2023    11  22   13 71142665      259.9323
## 6  2023-11-22 14:00:00      THYAO 2023    11  22   14 71142665      259.8585
## 7  2023-11-22 15:00:00      THYAO 2023    11  22   15 71142665      259.4808
## 8  2023-11-22 16:00:00      THYAO 2023    11  22   16 71142665      259.6106
## 9  2023-11-22 17:00:00      THYAO 2023    11  22   17 71142665      259.6953
## 10 2023-11-22 18:00:00      THYAO 2023    11  22   18 71142665      259.7274
##    price
## 1     NA
## 2     NA
## 3     NA
## 4     NA
## 5     NA
## 6     NA
## 7     NA
## 8     NA
## 9     NA
## 10    NA

Now, we have the data let’s do the predictions for all 30 stocks.

noftrees=500
depth=5
learning_rate=0.1
# Function to train random forest model and calculate performance metrics
train_and_evaluate_rf_model <- function(stock_data) {
  stock_name <- unique(stock_data$short_name)
  
  stock_data <- stock_data[, !colnames(stock_data) %in% c("short_name","timestamp")]

  validationIndex <- createDataPartition(stock_data$price, p=0.70, list=FALSE)
  
  gbm_train <- stock_data[validationIndex,] # 70% of data to training
  gbm_test <- stock_data[-validationIndex,]
  
  #ts_data <- ts(stock_data$price, frequency = 10)

  #arima_model <- auto.arima(ts_data)
  
  #predictions <- forecast(arima_model, h = nrow(rf_test))
  
  # Train GBM model
  boosting_model=gbm(price~., data=stock_data,distribution = 'gaussian', n.trees = noftrees,
                   interaction.depth = depth, n.minobsinnode = 10, shrinkage =learning_rate, cv.folds = 10)
  
  # Make predictions on test data
  predictions <- predict(boosting_model, newdata = gbm_test, type = "response")
  
  # Calculate performance metrics
  mae <- mae(predictions, gbm_test$price)
  mape <- mape(predictions, gbm_test$price)
  rmse <- rmse(predictions, gbm_test$price)
  wmape <- calculate_wmape(gbm_test$price, predictions)
  #mae <- mae(predictions$mean, rf_test$price)
  #mape <- mape(predictions$mean, rf_test$price)
  #rmse <- rmse(predictions$mean, rf_test$price)
  #wmape <- calculate_wmape(rf_test$price, predictions$mean)
  
  # Create a data frame to store performance metrics
  performance_df <- data.frame(StockName = stock_name, MAE = mae, RMSE = rmse, MAPE = mape, WMAPE = wmape)
  
  further_data <- new_data_df%>% filter(short_name == stock_name)
  # Predictions for further day
  further_predictions <- predict(boosting_model, newdata = further_data, type = "response")
  #further_predictions <- forecast(arima_model, h = nrow(further_data))

  assign(paste0("pred_", stock_name), further_predictions, envir = .GlobalEnv)
  
  return(performance_df)
}

# List to store performance metrics for each stock
performance_list <- list()

# List of stock data tables
stock_data_list <- list(thyao_data, akbnk_data, arclk_data, asels_data, bimas_data,
                        dohol_data, ekgyo_data, eregl_data, froto_data, gubrf_data,
                        garan_data, krdmd_data, kchol_data, kozal_data, kozaa_data,
                        pgsus_data, petkm_data, sahol_data, sasa_data, sise_data,
                        tahvl_data, tkfen_data, tuprs_data, ttkom_data, tcell_data,
                        halkb_data, isctr_data, vakbn_data, vestl_data, ykbnk_data)

# Train, evaluate, and store performance metrics for each stock
for (stock_data in stock_data_list) {
  performance_df <- train_and_evaluate_rf_model(stock_data)
  performance_list[[length(performance_list) + 1]] <- performance_df
}

# Combine performance metrics for all stocks into a single dataframe
all_performance_df <- do.call(rbind, performance_list)

all_performance_df

##    StockName        MAE       RMSE        MAPE     WMAPE
## 1      THYAO 0.41666947 0.81082198 0.010537753 0.8622404
## 2      AKBNK 0.08137366 0.13206826 0.009168904 0.9178386
## 3      ARCLK 0.39396395 0.68104401 0.009248987 0.8737846
## 4      ASELS 0.10763006 0.18730054 0.008912143 0.8996698
## 5      BIMAS 0.57860136 0.94768345 0.007499711 0.7402067
## 6      DOHOL 0.03442372 0.05943560 0.010993457 0.9845501
## 7      EKGYO 0.02655112 0.04632598 0.009378318 0.9294231
## 8      EREGL 0.13965642 0.22884386 0.008862340 0.8142383
## 9      FROTO 1.75369392 3.08642459 0.010191791 0.8485934
## 10     GUBRF 0.83944462 1.57674178 0.018067704 1.0648366
## 11     GARAN 0.12734486 0.22327770 0.009375542 0.9410761
## 12     KRDMD 0.08454021 0.14325373 0.011408820 1.0543903
## 13     KCHOL 0.29175883 0.49757475 0.008864704 0.8444175
## 14     KOZAL 0.07692361 0.15118136 0.010833130 1.0391935
## 15     KOZAA 0.21815909 0.37421116 0.011188007 1.0827450
## 16     PGSUS 1.58947450 2.97824774 0.012811443 0.9633349
## 17     PETKM 0.06531191 0.11116072 0.008759580 0.8859920
## 18     SAHOL 0.13075817 0.22411066 0.008324274 0.7983691
## 19      SASA 0.15119990 0.32418239 0.016288741 1.1172268
## 20      SISE 0.12042145 0.21318860 0.008943390 0.8191455
## 21     TAVHL 0.34852328 0.54769387 0.009522216 0.8895338
## 22     TKFEN 0.18086629 0.27914591 0.008788029 0.8939932
## 23     TUPRS 0.24444193 0.42469739 0.008731905 0.7994325
## 24     TTKOM 0.07981287 0.13116003 0.009523494 0.9375231
## 25     TCELL 0.15647807 0.24964013 0.008309426 0.8418590
## 26     HALKB 0.06391275 0.09981369 0.008620749 0.8857906
## 27     ISCTR 0.04678277 0.08529747 0.009711777 0.9664324
## 28     VAKBN 0.05411899 0.08561472 0.008958584 0.9239730
## 29     VESTL 0.23083211 0.38947893 0.010760648 0.9943066
## 30     YKBNK 0.04229955 0.07517906 0.010157310 0.9717961

ggplot(all_performance_df, aes(x = StockName)) +
  geom_bar(aes(y = MAE), fill = "blue", stat = "identity", position = "dodge") +
  geom_bar(aes(y = RMSE), fill = "green", stat = "identity", position = "dodge") +
  geom_bar(aes(y = MAPE), fill = "orange", stat = "identity", position = "dodge") +
  geom_bar(aes(y = WMAPE), fill = "red", stat = "identity", position = "dodge") +
  labs(title = "Performance Metrics by Stock",
       x = "Stock Name",
       y = "Value") +
  theme_minimal()

# Scatter plot for MAE and RMSE
ggplot(all_performance_df, aes(x = MAE, y = RMSE, color = StockName)) +
  geom_point() +
  labs(title = "Scatter Plot of MAE vs RMSE",
       x = "Mean Absolute Error (MAE)",
       y = "Root Mean Squared Error (RMSE)") +
  theme_minimal()

Stock Prices Submission Generator

Since we are going to submit our predictions through Google Forms as a pre-defined format. We wrote a function that creates this predictions submission format. It iteratively adds prediction (e.g. pred_THYAO, pred_SISE etc.) to a key as stock name in a dictionary.

stock_submission_generator <- function(stock_names) {
  # Create a list to store the predictions
  prediction_list <- list()
  
  # Iterate over each stock name
  for (stock_name in stock_names) {
    # Clean up the stock name to make it a valid R variable name
    cleaned_stock_name <- gsub("[^A-Za-z0-9_]", "_", stock_name)
    
    # Get the predicted values directly (no need for dynamic variable names)
    prediction_values <- get(paste0("pred_", cleaned_stock_name))
    
    # Store the results in the list
    prediction_list[[cleaned_stock_name]] <- prediction_values
  }
  # Combine the results into a named list
    result <- paste0("{", paste(lapply(names(prediction_list), function(stock_name) {
      paste0("'", stock_name, "': [", paste(prediction_list[[stock_name]], collapse = ", "), "]")
    }), collapse = ","), "}")
  
  return(result)
}

# Example usage
stock_names <- c("THYAO", "AKBNK", "ARCLK", "ASELS", "BIMAS", "DOHOL", "EKGYO", "EREGL", "FROTO", "GUBRF",
                 "GARAN", "KRDMD", "KCHOL", "KOZAL", "KOZAA", "PGSUS", "PETKM", "SAHOL", "SASA", "SISE",
                 "TAVHL", "TKFEN", "TUPRS", "TTKOM", "TCELL", "HALKB", "ISCTR", "VAKBN", "VESTL", "YKBNK")

# Call the function with the stock names and their respective predictions
submission_result <- stock_submission_generator(stock_names)

cat(submission_result)

## {'THYAO': [259.678273415031, 259.488180755706, 262.870441901668, 264.575878798915, 264.928903323283, 263.854646405149, 262.843440688276, 262.93481263193, 263.280647923134, 263.326349039752],'AKBNK': [30.7488662116627, 30.6865904971919, 30.6543959667005, 30.539026700447, 30.7834293207161, 30.7487436661858, 30.4628626078624, 30.4763084611589, 30.3906350771, 30.3876187209088],'ARCLK': [148.42956749188, 148.904498176351, 149.237481143448, 149.325499844081, 148.822245211408, 148.677481699621, 147.59928510894, 147.493642461175, 146.803227009647, 146.630261067425],'ASELS': [44.4421644435841, 44.4758067658741, 44.5470324433092, 44.4802434325592, 44.5145429357864, 44.5626216287134, 44.6675513715604, 44.6945284241698, 44.7393691194361, 44.8023816606475],'BIMAS': [297.462952177863, 298.436158960658, 301.213758714903, 302.126539418559, 303.835168271841, 304.900077047282, 304.308511211628, 304.876119961817, 305.002624147136, 305.002624147136],'DOHOL': [13.6736249567136, 13.6822175384218, 13.7019084840272, 13.6899924993666, 13.6713149226958, 13.6509629172961, 13.6488663142748, 13.7028655126324, 13.6997877393759, 13.7093936219274],'EKGYO': [7.70686192514752, 7.7249728045728, 7.72432586379302, 7.70808691466281, 7.70097680283151, 7.67482832011071, 7.65647475724637, 7.63342174538909, 7.62768202262186, 7.62933857424417],'EREGL': [40.7114524054011, 40.9521147856662, 40.9502447892473, 40.9095586076127, 40.8748902880372, 40.8275925402601, 40.7907905406097, 40.7946364986888, 40.7677466071393, 40.7742824326827],'FROTO': [864.726999045236, 867.903508478008, 869.53949905273, 869.29845669823, 867.682615902287, 866.794276618594, 864.407182528025, 867.059894193871, 865.814762723779, 865.586728665435],'GUBRF': [377.61937252875, 375.301927618272, 374.412949821756, 374.277121603163, 375.589118204399, 376.868324559247, 376.631428603911, 376.434489838284, 375.970998250477, 375.70034328655],'GARAN': [49.247985848862, 49.3631813934211, 49.4522897124222, 49.3421178946062, 49.3173476374261, 49.2875353281418, 49.2748924870698, 49.3634975554617, 49.465448480069, 49.508729677727],'KRDMD': [25.0040465576221, 25.1130929130406, 25.1070580340409, 25.0794698485084, 25.051184854409, 25.01426598558, 24.9910636924458, 24.9756322444879, 24.985484330764, 25.0035025458213],'KCHOL': [142.630104942437, 142.763542880071, 142.959672068049, 143.585174962454, 143.547306163818, 143.67368486791, 142.527112631001, 142.350140233165, 142.028231009, 141.650921230105],'KOZAL': [23.3601767478346, 23.329664137582, 23.3415675220427, 23.3055368691179, 23.3059745480377, 23.1824384869305, 23.1195886650377, 23.1105491407486, 23.0719007966935, 23.0719007966935],'KOZAA': [57.5323460638849, 57.4415318967575, 56.6919598601638, 57.0734397473845, 56.9491658029345, 56.9775620310227, 56.5421008698111, 56.1076664521586, 55.4256299166737, 55.4586487002782],'PGSUS': [741.768993415986, 740.722006679985, 739.184357710692, 739.24259818296, 737.948238704466, 738.272142315949, 736.800397180212, 738.988978984747, 739.637224967399, 739.223316185482],'PETKM': [22.5060575546983, 22.4630456908606, 22.4420915790538, 22.1918589394313, 22.1732111799991, 22.3886778981182, 22.3680260231842, 22.4566025239067, 22.4535935297212, 22.4482486376441],'SAHOL': [59.7838994560353, 59.8302846634545, 60.118679680521, 60.2268200003305, 60.2918521541968, 60.2165632502911, 60.1771180067082, 60.1637058359442, 60.1023126467572, 60.2082496431509],'SASA': [45.6422005746832, 45.6991144364409, 45.6832712644442, 45.660340467497, 45.660340467497, 45.591926008355, 45.3753743526816, 45.3753743526816, 45.2674663919705, 45.2674663919705],'SISE': [50.6454575072272, 50.7373348779948, 50.8135046025052, 50.7275202594185, 50.6825241583468, 50.657595864405, 50.6939335389986, 50.7060544218619, 50.6697054591102, 50.7015412637805],'TAVHL': [126.07689615251, 126.043956590835, 125.849141670105, 125.696886120224, 125.33213125202, 125.350036622323, 125.120910207309, 124.834802566218, 125.04945030825, 125.260213260965],'TKFEN': [46.0763239075224, 46.0293046530382, 46.1728452793203, 46.2025532900501, 46.2766168067448, 46.3803884624692, 46.4056328050447, 46.2925182428383, 46.1813799621897, 46.1993263429708],'TUPRS': [156.098424106778, 156.247377827694, 156.262271902691, 155.771414452609, 155.571838532015, 154.930006363034, 154.628723346755, 154.510708912732, 154.021628744227, 154.074929950634],'TTKOM': [22.3931698843898, 22.3805528990429, 22.3635538613355, 22.3474631154122, 22.3019110575896, 22.2684444574148, 22.2479566447482, 22.2677927667868, 22.2607729282805, 22.3002616457167],'TCELL': [56.9190228694598, 56.9980248863929, 57.0173768445357, 56.982589641217, 56.8168684766495, 56.5892589934113, 56.3706496624305, 56.289961235087, 55.9721628512774, 56.0227791591829],'HALKB': [13.1798344122416, 13.1815902545908, 13.1811449054022, 13.1653038058198, 13.1638255987007, 13.1390001752411, 13.1396112986481, 13.1285456685077, 13.1018692293052, 13.1018692293052],'ISCTR': [20.6047874413227, 20.6260827012082, 20.6130046479076, 20.6476735570818, 20.6532750154644, 20.6554618493375, 20.6544189281589, 20.6463566099548, 20.6365691036402, 20.6365691036402],'VAKBN': [15.1140002196792, 15.0986012053951, 15.0996388053073, 15.1066729967723, 15.1101395681491, 15.0839499890343, 15.0671565986052, 15.0715368771805, 15.0216823920686, 15.031657001491],'VESTL': [63.3875065033263, 63.7297887697724, 63.7614448809427, 63.5000986786597, 63.2610077519992, 62.9444828063781, 62.6825178878897, 62.5501564057271, 62.1785786910445, 62.3097569111928],'YKBNK': [17.0965708995891, 17.1153016007646, 17.1153016007646, 17.0953543802834, 17.08794997852, 17.0905240719484, 17.0898712804158, 17.1087461139085, 17.0744953415476, 17.0856990295545]}

Calculating Next Day’s Expected Volume

Since Yahoo Finance Library enables volume information of stocks until day t and we need to predict next days’ stock prices. We have written a function that returns the excepted volume of a specific stock of the next day. It basically, aggregates 50 previous days and returns the calculate of their average volume.

calculate_next_day_expected_volume <- function(stock_name){
  
  test <- final_data %>% filter(short_name == stock_name)
  
  expected_nextday_volume <- round(mean(tail(test$volume, 50)))
  
  day <- tail(test$timestamp,1)
  
  timestamp <- ymd_hms(day, tz = "UTC")
  
  # Increment the date by 1 day and floor it to remove hours
  next_day <- floor_date(timestamp + days(1), unit = "days")
  
  cat("Volume for the stock",stock_name," for the day ",format(next_day, "%Y-%m-%d"),"is :", expected_nextday_volume)
  
}

calculate_next_day_expected_volume("THYAO")

## Volume for the stock THYAO  for the day  2023-11-22 is : 47201751

calculate_next_day_expected_volume("AKBNK")

## Volume for the stock AKBNK  for the day  2023-11-22 is : 68234660

calculate_next_day_expected_volume("ARCLK")

## Volume for the stock ARCLK  for the day  2023-11-22 is : 3345512

calculate_next_day_expected_volume("FROTO")

## Volume for the stock FROTO  for the day  2023-11-22 is : 791529

Predicting Stock Prices with Unseen Data

In this part, we are going to see that how well our selected model predicts the unseen next day’s stock price.

new_data <- data.frame(
  timestamp = seq(as.POSIXct("2023-12-12 09:00:00"), by = "hour", length.out = 5),
  short_name = rep("THYAO", 5),
  year = c(2018, 2018, 2018, 2018, 2018),
  month = c(1, 1, 1, 1, 1),
  day = c(2, 2, 2, 2, 2),
  hour = c(9, 10, 11, 12, 13),
  volume = rep(36869621 , 5),
  forecast_ma_6 = c(260.2, 260.8, 265.2, 262.5, 265),
  price = c(261.5, 262.4, 263.7, 264, 263)
)

new_data$year[1:5] <- 2023
new_data$month[1:5] <- 12
new_data$day[1:5] <- 12
new_data$hour[1:5] <- 9:13


head(new_data)

##             timestamp short_name year month day hour   volume forecast_ma_6
## 1 2023-12-12 09:00:00      THYAO 2023    12  12    9 36869621         260.2
## 2 2023-12-12 10:00:00      THYAO 2023    12  12   10 36869621         260.8
## 3 2023-12-12 11:00:00      THYAO 2023    12  12   11 36869621         265.2
## 4 2023-12-12 12:00:00      THYAO 2023    12  12   12 36869621         262.5
## 5 2023-12-12 13:00:00      THYAO 2023    12  12   13 36869621         265.0
##   price
## 1 261.5
## 2 262.4
## 3 263.7
## 4 264.0
## 5 263.0

With Decision Tree:

test_pred <- predict(best_tree_model, new_data[,-9])
test_pred

##      1      2      3      4      5 
## 259.04 259.04 259.04 259.04 259.04

pred_THYAO <- predict(boosting_model, new_data[,-9])
mae <- mean(abs(test_pred - new_data$price))
print(paste("Mean Absolute Error (MAE):", mae))

## [1] "Mean Absolute Error (MAE): 3.87999999999997"

# Calculate Root Mean Squared Error (RMSE)
rmse <- round(sqrt(mean((test_pred - new_data$price)^2)),3)
print(paste("Root Mean Squared Error (RMSE):", rmse))

## [1] "Root Mean Squared Error (RMSE): 3.983"

# Calculate Mean Absolute Percentage Error (MAPE)
mape <- mean(abs((test_pred - new_data$price) / new_data$price) * 100, na.rm = TRUE)
cat("Mean Absolute Percentage Error (MAPE):", mape, "\n")

## Mean Absolute Percentage Error (MAPE): 1.474573

# Calculate Weighted Mean Absolute Percentage Error (WMAPE)
wmape <- calculate_wmape(test_pred, new_data$price)
print(paste("Weighted Mean Absolute Percentage Error (WMAPE):", wmape))

## [1] "Weighted Mean Absolute Percentage Error (WMAPE): 1.49783817171092"

With ARIMA:

forecast_values <- forecast(arima_model, h = nrow(new_data))  
forecast_values$mean

## Time Series:
## Start = c(1025, 8) 
## End = c(1026, 2) 
## Frequency = 10 
## [1] 270.9827 273.1754 275.6821 278.4079 281.3058

mae <- mean(abs(forecast_values$mean - new_data$price))
mse <- mean((forecast_values$mean - new_data$price)^2)
rmse <- sqrt(mse)
mape <- mean(abs((forecast_values$mean - new_data$price) / new_data$price) * 100, na.rm = TRUE)
wmape <- calculate_wmape(forecast_values$mean, new_data$price)
cat("MAE:", mae, "\n")

## MAE: 12.99077

cat("MSE:", mse, "\n")

## MSE: 178.458

cat("RMSE:", rmse, "\n")

## RMSE: 13.35882

cat("MAPE:", mape, "\n")

## MAPE: 4.938897

cat("WMAPE:",wmape, "\n")

## WMAPE: 4.708322